import pandas as pd
import numpy as np
import geopandas as gpd
import contextily as cx
import matplotlib.pyplot as plt
from shapely.geometry import Point, MultiPoint, MultiPolygon, Polygon, LineString, MultiLineString
import networkx as nx
import momepy as mp
import tqdm
import libpysal

# OS ROADS # 
OS_rds =  gpd.read_file("C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/ordinance_survey_roads_BBs.gpkg")
OS_rds = OS_rds[['city', 'geometry']]
OS_rds = gpd.GeoDataFrame(OS_rds)
OS_rds = OS_rds.explode()

lrds = OS_rds.loc[OS_rds['city']=='Leicester']

lrds  = mp.remove_false_nodes(lrds)
lrds.geometry = mp.close_gaps(lrds, 2)
lrds['city'] = 'Leicester'
lrds = lrds[['geometry']]
lrds['nID'] = mp.unique_id(lrds)

# # # CITY BOUNDARIES #
boundaries = gpd.read_file("C:/Users/patri/Documents/CDRC Bikeshare Station Data/spatial_data_pre_process/combined_boundaries_updated_poole_christchurch_clipped.gpkg")
boundaries = boundaries.rename(columns={'city':'boundary_city'})
boundaries['boundary_city'] = boundaries['boundary_city'].replace({'Glasgow City': 'Glasgow', 'City of Leicester': 'Leicester'})
boundaries = boundaries.loc[boundaries['boundary_city']=='Leicester']

# BUILDINGS  # 
lbuildings = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/leicester_buildings_OS.gpkg')

# BUILT UP AREAS # 
built_up_areas = gpd.read_file("C:/Users/patri/Downloads/OS_Open_Built_Up_Areas_GeoPackage/OS_Open_Built_Up_Areas.gpkg")
built_up_areas.to_crs(27700, inplace=True)
built_up_areas = built_up_areas[['geometry']]

def building_id(buildings):
    buildings = buildings.to_crs(27700)
    buildings['uID'] = mp.unique_id(buildings)
    buildings = buildings[['uID', 'geometry']]

    return buildings

lbuildings = building_id(lbuildings)

# TESSELATION FUNCTION # 

def tesselate_func(buildings, streets_gdf, letter):

    limit = mp.buffered_limit(buildings, 100)

    tessellation = mp.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
    tessellation = tessellation.tessellation
    print("tesseltation DONE")
    tessellation = tessellation.overlay(built_up_areas, how='intersection')
    
    print("tesseltation INTERSECTION DONE")

    #link streets
    buildings = buildings[['uID', 'geometry']]
    buildings = buildings.sjoin_nearest(streets_gdf, max_distance=1000)
    buildings = buildings.drop_duplicates("uID").drop(columns="index_right")
    tessellation = tessellation.merge(buildings[['uID', 'nID']], on='uID', how='left')
    tessellation = tessellation.drop_duplicates("uID")
    
    # Dimensions
    buildings["area"] = buildings.area
    tessellation["area"] = tessellation.area
    streets_gdf["length"] = streets_gdf.length
    # Shape
    buildings['eri'] = mp.EquivalentRectangularIndex(buildings).series
    print('EquivalentRectangularIndex DONE')
    buildings['elongation'] = mp.Elongation(buildings).series
    print('elongation DONE')
    tessellation['convexity'] = mp.Convexity(tessellation).series
    print('convexity DONE')
    streets_gdf["linearity"] = mp.Linearity(streets_gdf).series
    print('linearity DONE')

    # Spatial distribution
    buildings["shared_walls"] = mp.SharedWallsRatio(buildings).series
    print('shared_walls DONE')

    # Generate spatial weights matrix using libpysal.

    queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)
    tessellation["neighbors"] = mp.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
    tessellation["covered_area"] = mp.CoveredArea(tessellation, queen_1, "uID", verbose=False).series

    buildings["neighbor_distance"] = mp.NeighborDistance(buildings, queen_1, "uID", verbose=False).series

    queen_3 = mp.sw_high(k=3, weights=queen_1)
    buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)

    # buildings['interbuilding_distance'] = mp.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
    buildings['adjacency'] = mp.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
    print('adjacency DONE') 
    buildings['squarenes'] = mp.Squareness(buildings).series
    print('squarenes DONE') 
    buildings['squarecompactness'] = mp.SquareCompactness(buildings).series
    print('squarecompactness DONE') 
    buildings['corners'] = mp.Corners(buildings).series
    print('corners DONE') 

    profile = mp.StreetProfile(streets_gdf, buildings)
    streets_gdf["width"] = profile.w
    print('width DONE') 
    streets_gdf["width_deviation"] = profile.wd
    print('width_deviation DONE') 
    streets_gdf["openness"] = profile.o
    print('openness DONE')

    # Intensity
    tessellation['area_ratio'] = mp.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
    print('area_ratio DONE')

    # Connectivity
    graph = mp.gdf_to_nx(streets_gdf)

    graph = mp.node_degree(graph)
    print('node degree DONE')
    graph = mp.closeness_centrality(graph, radius=2000, distance="mm_len")
    print('closeness DONE')
    graph = mp.meshedness(graph, radius=2000, distance="mm_len")
    print('meshedness DONE')
    graph = mp.gamma(graph, radius=2000, name='gamma', distance="mm_len", verbose=True)
    print('gamma DONE')
    graph = mp.edge_node_ratio(graph, radius=500, name='edge_node_ratio', distance="mm_len", verbose=True)
    print('edge node ratio DONE')
    nodes, streets = mp.nx_to_gdf(graph)

    streets['street_orientation'] = mp.NeighboringStreetOrientationDeviation(streets).series
    print('street_orientation ratio DONE')

    buildings["nodeID"] = mp.get_node_id(buildings, nodes, streets, "nodeID", "nID")
    # Link all data together (to tessellation cells or buildings).

    # tessellation
    streets = streets[['nID', 'street_orientation','geometry']]
    merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
    merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
    merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')

    return buildings, streets, nodes, tessellation, merged

lbuildings, lstreets, lnodes, ltessellation, lmerged = tesselate_func(lbuildings, lrds, "L")

merged = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/L_merged_tessellation_100m.gpkg', driver='GPKG')
buildings = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/L_buildings_tesselate_stats_100m.gpkg', driver='GPKG')
streets = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/L_streets_gdf_tesselate_stats_100m.gpkg', driver='GPKG')
tessellation = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/L_tesselation_100m.gpkg', driver='GPKG')
nodes = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/L_nodes_tesselation_100m.gpkg', driver='GPKG')

nodes.columns

Index(['degree', 'cds_len', 'closeness2000', 'meshedness', 'gamma2000',
       'edge_node_ratio', 'nodeID', 'closeness_global', 'geometry'],
      dtype='object')

nodes

fig, axs = plt.subplots(2, 4, figsize=(32, 16))

# Define a function to style each subplot
def style_subplot(ax, title):
    ax.set_title(title, fontsize=16)
    ax.axis('off')  # This removes all axis lines, ticks, and labels

# Plot each subplot with reduced point size and custom styling
variables = [
    "closeness2000", "closeness_global", "meshedness", "degree",
    "cds_len", "edge_node_ratio", "gamma2000"
]

for i, var in enumerate(variables):
    row = i // 4
    col = i % 4
    if var:  # Only plot if there's a variable name
        nodes.plot(var, ax=axs[row, col], scheme='natural_breaks', legend=True, markersize=5)
        style_subplot(axs[row, col], var.replace("_", " ").title())
    else:
        axs[row, col].axis('off')  # Turn off empty subplots

# Adjust layout to prevent overlap
plt.tight_layout()

# Show the plot
plt.show()

streets = gpd.clip(streets, boundaries)

fig, axs = plt.subplots(2, 3, figsize=(24, 16))

# Define a function to style each subplot
def style_subplot(ax, title):
    ax.set_title(title, fontsize=16)
    ax.axis('off')  # This removes all axis lines, ticks, and labels

# Plot each subplot with reduced point size and custom styling
variables = [
    "linearity", "length", "width", "width_deviation",
    "continuity", "openness"
]

for i, var in enumerate(variables):
    row = i // 3
    col = i % 3
    if var:  # Only plot if there's a variable name
        streets.plot(var, ax=axs[row, col], scheme='natural_breaks', legend=True, markersize=5)
        style_subplot(axs[row, col], var.replace("_", " ").title())
    else:
        axs[row, col].axis('off')  # Turn off empty subplots

# Adjust layout to prevent overlap
plt.tight_layout()

# Show the plot
plt.show()

buildings.columns

Index(['uID', 'nID', 'area', 'eri', 'elongation', 'neighbor_distance',
       'shared_walls', 'swr', 'interbuilding_distance', 'adjacency',
       'squarenes', 'squarecompactness', 'corners', 'comp',
       'building_orientation', 'cell_alignment', 'nodeID', 'geometry'],
      dtype='object')

fig, axs = plt.subplots(2, 3, figsize=(24, 16))

# Define a function to style each subplot
def style_subplot(ax, title):
    ax.set_title(title, fontsize=16)
    ax.axis('off')  # This removes all axis lines, ticks, and labels

# Plot each subplot with reduced point size and custom styling
variables = [
    "eri", "elongation", "neighbor_distance",  'building_orientation',
    "shared_walls", "squarecompactness"
]

for i, var in enumerate(variables):
    row = i // 3
    col = i % 3
    if var:  # Only plot if there's a variable name
        buildings.plot(var, ax=axs[row, col], scheme='natural_breaks', legend=True)
        style_subplot(axs[row, col], var.replace("_", " ").title())
    else:
        axs[row, col].axis('off')  # Turn off empty subplots

# Adjust layout to prevent overlap
plt.tight_layout()

# Show the plot
plt.show()

from shapely.geometry import box
# BBS
city_bbs = gpd.read_file("C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/bounding_boxes_stations_boundaries.gpkg")
city_bbs.to_crs(4326, inplace=True)
lbb = city_bbs.loc[city_bbs['city']=='Leicester']
lbb.to_crs(27700, inplace=True)
lbb['centroid'] = lbb.geometry.centroid
lbb.set_geometry('centroid', inplace=True)

c:\Users\patri\Documents\spopt_venv\lib\site-packages\geopandas\geodataframe.py:1528: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
c:\Users\patri\Documents\spopt_venv\lib\site-packages\geopandas\geodataframe.py:1528: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)

limit = mp.buffered_limit(buildings, 100)
l_tess = mp.Tessellation(buildings, 'uID', limit, segment=1) # create tesselations from all buildings
l_tess = l_tess.tessellation

# Create the bounding box
def create_bbox(point, distance):
    minx, miny, maxx, maxy = point.buffer(distance).bounds
    return box(minx, miny, maxx, maxy)

# Apply the function to create a 300m x 300m bounding box
distance = 150  # 150 meters in each direction to create a 300m x 300m box
lbb['bbox'] = lbb.geometry.apply(lambda point: create_bbox(point, distance))

# Create a new GeoDataFrame with the bounding boxes
bbox_gdf = gpd.GeoDataFrame(geometry=lbb['bbox'], crs=lbb.crs)

# Filter buildings and roads to the bounding box
buildings_in_bbox = gpd.overlay(buildings, bbox_gdf)
buildings_in_bbox_buff = buildings_in_bbox.buffer(20)
buildings_in_bbox_buff = gpd.GeoDataFrame(buildings_in_bbox_buff)
buildings_in_bbox_buff.rename(columns={0:'geometry'}, inplace=True)
buildings_in_bbox_buff.set_geometry('geometry', inplace=True)
buildings_in_bbox_buff = buildings_in_bbox_buff.dissolve()
roads_in_bbox = gpd.overlay(streets, bbox_gdf)
l_tess_in_bbox = gpd.overlay(l_tess, buildings_in_bbox_buff)

Inward offset...
Generating input point array...
Generating Voronoi diagram...
Generating GeoDataFrame...
Dissolving Voronoi polygons...

C:\Users\patri\AppData\Local\Temp\ipykernel_11772\2741401941.py:2: UserWarning: Tessellation does not fully match buildings. 10 element(s) collapsed during generation - unique_id: {25986, 48455, 26379, 26444, 51162, 48348, 48483, 26349, 48365, 2230}.
  l_tess = mp.Tessellation(buildings, 'uID', limit, segment=1)
C:\Users\patri\AppData\Local\Temp\ipykernel_11772\2741401941.py:2: UserWarning: Tessellation contains MultiPolygon elements. Initial objects should  be edited. `unique_id` of affected elements: [25863, 3296].
  l_tess = mp.Tessellation(buildings, 'uID', limit, segment=1)
c:\Users\patri\Documents\spopt_venv\lib\site-packages\geopandas\geodataframe.py:1528: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[38], line 24
     22 buildings_in_bbox_buff.set_geometry('geometry', inplace=True)
     23 buildings_in_bbox_buff = buildings_in_bbox_buff.dissolve()
---> 24 roads_in_bbox = gpd.overlay(lrds, bbox_gdf)
     25 l_tess_in_bbox = gpd.overlay(ltessellation, buildings_in_bbox_buff)

NameError: name 'lrds' is not defined

# Plot the results
fig, ax = plt.subplots(2, 2, figsize=(10, 10))
#1 
roads_in_bbox.plot(ax=ax[0, 0], color='grey', linewidth=1)
#2
buildings_in_bbox.plot(ax=ax[0, 1], color='gray', alpha=0.5)
#3
buildings_in_bbox.plot(ax=ax[1, 0], color='gray', alpha=0.5)
roads_in_bbox.plot(ax=ax[1, 0], color='grey', linewidth=1)
#4
buildings_in_bbox.plot(ax=ax[1,  1], color='gray', alpha=0.5)
roads_in_bbox.plot(ax=ax[1, 1], color='grey', linewidth=1, alpha=0.7)
l_tess_in_bbox.plot(ax=ax[1,1], linewidth=1.2, edgecolor='blue', alpha=0.5)
buildings_in_bbox_buff.plot(ax=ax[1,1], alpha=0.1, color='gray', edgecolor='grey')

ax[0,0].set_title("A) Road segments", fontsize=12)
ax[0,1].set_title("B) Building polygons", fontsize=12)
ax[1,0].set_title("C) Roads and building polygons joined", fontsize=12)
ax[1,1].set_title("D) Tessellated units, bounded by a distance-band", fontsize=12)
ax[0,0].axis('off') 
ax[0,1].axis('off')  
ax[1,0].axis('off')  
ax[1,1].axis('off')

(458826.215508731, 459200.215508731, 304404.94278221956, 304778.94278221956)

fig, axs = plt.subplots(2, 3, figsize=(24, 16))

# Define a function to style each subplot
def style_subplot(ax, title):
    ax.set_title(title, fontsize=16)
    ax.axis('off')  # This removes all axis lines, ticks, and labels

# Plot each subplot with reduced point size and custom styling
variables = [
    "closeness2000", "cds_len", "neighbor_distance",  'linearity',
    "width_deviation", "squarecompactness"
]

for i, var in enumerate(variables):
    row = i // 3
    col = i % 3
    if var:  # Only plot if there's a variable name
        merged.plot(var, ax=axs[row, col], scheme='natural_breaks', legend=True)
        style_subplot(axs[row, col], var.replace("_", " ").title())
    else:
        axs[row, col].axis('off')  # Turn off empty subplots

# Adjust layout to prevent overlap
plt.tight_layout()

# Show the plot
plt.show()

from libpysal.weights import DistanceBand
om500 = DistanceBand.from_dataframe(merged, 500, geom_col='centroid', ids='xID')

def create_iqr_columns(df, columns, weight_matrix, id_col, rng=(25, 75)):
    
    columns_to_process = [col for col in df.columns if col != id_col]
    
    for column in columns_to_process:
        new_column_name = f"{column}_IQR_db500"
        print(f"Processing column: {column}")
        print(f"Weight matrix: {weight_matrix}")
        
        try:
            df[new_column_name] = mp.Range(df, column, weight_matrix, id_col, rng=rng).series
        except Exception as e:
            print(f"Error processing column {column}: {e}")

IQR_db500_df = create_iqr_columns(merged, cluster_cols, om500, 'xID') # I have ommited the cluster columns originally included for brevity.

cluster_data = pd.read_csv('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/cluster_data_distance_bands_100m_buffer_2.csv')

# GEOM
merged_all_geom = gpd.read_file('C:/Users/patri/Documents/CDRC Bikeshare Station Data/station_optimisation/data_pre_process/tessellation/merged_all_geom_100m_buffer_2.gpkg', driver='GPKG')
IQR_db500_df = pd.concat( [cluster_data,merged_all_geom[['city','geometry']]], axis=1)
IQR_db500_df = gpd.GeoDataFrame(IQR_db500_df, geometry='geometry')

IQR_db500_df = pd.concat( [cluster_data,merged_all_geom[['geometry']]], axis=1)
IQR_db500_df = gpd.GeoDataFrame(IQR_db500_df, geometry='geometry')
IQR_db500_df = IQR_db500_df.loc[IQR_db500_df['city']=='Leicester']

fig, axs = plt.subplots(2, 3, figsize=(24, 16))

# Define a function to style each subplot
def style_subplot(ax, title):
    ax.set_title(title, fontsize=16)
    ax.axis('off')  # This removes all axis lines, ticks, and labels

# Plot each subplot with reduced point size and custom styling
variables = [
    "closeness2000_IQR_db500", "cds_len_IQR_db500", "neighbor_distance_IQR_db500",  'linearity_IQR_db500',
    "width_deviation_IQR_db500", "squarecompactness_IQR_db500"
]

for i, var in enumerate(variables):
    row = i // 3
    col = i % 3
    if var:  # Only plot if there's a variable name
        IQR_db500_df.plot(var, ax=axs[row, col], scheme='natural_breaks', legend=True)
        style_subplot(axs[row, col], var.replace("_", " ").title())
    else:
        axs[row, col].axis('off')  # Turn off empty subplots

# Adjust layout to prevent overlap
plt.tight_layout()

# Show the plot
plt.show()

	degree	cds_len	closeness2000	meshedness	gamma2000	edge_node_ratio	nodeID	closeness_global	geometry
0	3	15129.056047	0.000035	0.111384	0.408074	1.000000	0	0.000148	POINT (457073.270 300462.100)
1	3	15733.159494	0.000035	0.110479	0.407468	1.020000	1	0.000147	POINT (457070.000 300432.000)
2	3	24310.011935	0.000027	0.108954	0.406609	1.071429	2	0.000151	POINT (463239.000 303470.000)
3	3	22144.062804	0.000026	0.111880	0.408602	1.117647	3	0.000146	POINT (463314.000 303417.000)
4	4	17293.273810	0.000039	0.146743	0.431569	1.291667	4	0.000147	POINT (458455.000 300635.000)
...	...	...	...	...	...	...	...	...	...
13406	1	11783.593790	0.000016	0.101806	0.402186	0.800000	13406	0.000183	POINT (464448.000 305431.000)
13407	3	9593.188961	0.000021	0.090656	0.394615	1.074074	13407	0.000001	POINT (464808.000 304694.000)
13408	1	4625.099820	0.000002	0.044776	0.372549	0.875000	13408	0.000001	POINT (464906.280 298069.010)
13409	1	7022.881817	0.000016	0.111298	0.408530	0.923077	13409	0.000000	POINT (464927.010 305723.750)
13410	1	6282.603608	0.000007	0.061224	0.376694	0.933333	13410	0.000000	POINT (464960.110 303724.480)

Measuring urban morphology - a guide¶

Why?¶

Initialising the libraries, data and functions¶

Visualising the data¶

Checking the road network¶

Checking the road network measures from road segments¶

Checking the building footprints¶

Checking the tessellated units¶

Drawing meaningful insights from the data¶

Why have we done this?¶